This report runs base, AR1, and AR2 models (variable ~ time) for monthly SEV meteorological data using the nlme R package. Mean monthly air temperature, minimum monthly air temperature, maximum monthly air temperature, and total monthly precipitation are modeled. More sophisticated models will likely be required after this preliminary analysis. The main goal is to determine whether certain variables are changing over time for different months of the year for the SEV meteorological stations.

library(tidyverse)
library(lubridate)
library(nlme)
library(MuMIn)
library(gridExtra)
library(kableExtra)

# Load yearly data ---------------------------------------------

path_to_files <- "../data/processed_data/"

# load and only keep vars of interest
m <- read_csv(paste0(path_to_files, "met_monthly_gap_filled.csv")) %>% 
  mutate(sta = as.factor(sta)) %>% 
  select(sta:ppt)
str(m)
## tibble [1,320 × 8] (S3: tbl_df/tbl/data.frame)
##  $ sta   : Factor w/ 4 levels "40","42","49",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ month : num [1:1320] 1 1 1 1 1 1 1 1 1 1 ...
##  $ year  : num [1:1320] 1988 1989 1990 1991 1992 ...
##  $ date  : Date[1:1320], format: "1988-01-01" "1989-01-01" ...
##  $ airt  : num [1:1320] 2.291 2.291 2.291 0.615 -1.32 ...
##  $ maxair: num [1:1320] 11.9 11.9 11.9 12.6 15.9 ...
##  $ minair: num [1:1320] -7.2 -7.2 -7.2 -13.8 -15.1 ...
##  $ ppt   : num [1:1320] 0 0 3.17 7.2 16.43 ...
summary(m)
##  sta          month            year           date                 airt       
##  40:408   Min.   : 1.00   Min.   :1988   Min.   :1988-01-01   Min.   :-1.320  
##  42:396   1st Qu.: 3.75   1st Qu.:2001   1st Qu.:2001-02-22   1st Qu.: 6.362  
##  49:276   Median : 6.50   Median :2008   Median :2008-03-16   Median :13.724  
##  50:240   Mean   : 6.50   Mean   :2007   Mean   :2007-07-13   Mean   :13.889  
##           3rd Qu.: 9.25   3rd Qu.:2015   3rd Qu.:2015-02-01   3rd Qu.:21.895  
##           Max.   :12.00   Max.   :2021   Max.   :2021-12-01   Max.   :29.569  
##      maxair           minair             ppt         
##  Min.   : 8.985   Min.   :-29.930   Min.   :  0.000  
##  1st Qu.:21.475   1st Qu.: -9.523   1st Qu.:  4.104  
##  Median :28.725   Median : -2.680   Median : 13.524  
##  Mean   :28.059   Mean   : -1.336   Mean   : 21.875  
##  3rd Qu.:34.655   3rd Qu.:  7.400   3rd Qu.: 31.050  
##  Max.   :48.720   Max.   : 17.910   Max.   :183.147

Prepare data for modeling -

Split data into individual met stations and add a time variable for nlme modeling.

# function to select station and variable and add 'time' var for a specific month of the year
prepare_data_for_nlme <- function(data, sta, month) {
  # data = dataframe to use
  # sta = met station (in quotes)
  # month = month of year (in quotes)
  data %>% 
    filter(sta == {{ sta }} & month == {{ month }}) %>% 
    arrange(date) %>% 
    mutate(time = as.numeric(as.factor(year)))
}

Set up functions -

This code produces several functions that will be useful for producing cleaner code.

# function for plotting preliminary monthly graph -----------------------------------------------------------

plot_monthly_prelim <- function(data, var, title, y_axis_label) {
  # data = dataset
  # var = variable (not in quotes)
  # title = title for graph (in quotes)
  # y_axis_label = label for y-axis (in quotes)
  ggplot(data, aes(x = date, y = {{ var }}, col = sta)) +
    geom_line(color = "burlywood") +
    geom_smooth(method = "lm", color = "black", size = 0.6) +
    facet_wrap(~ sta) +
    labs(title = title,
         x = "Month",
         y = y_axis_label)
}



by_month_lm <- function(data, var, station, title, y_axis_label) {
  # data = dataset
  # var = variable to plot (not in quotes)
  # station = met station (in quotes)
  # title = graph title
  # y_axis_label = label for the graph's y-axis
  m %>% 
    filter(sta == station) %>% 
    ggplot(., aes(x = year, y = {{ var }})) +
    geom_line() +
    geom_smooth(method = "lm") +
    facet_wrap(~ month, scales = "free_y") +
    ggtitle({{ title }})
}


# function for plotting graphs with best nlme model results ------------------------------------------------

plot_monthly_results <- function(data, var, coeffs, title, y_axis_label) {
  # data = dataset
  # var = variable (not in quotes)
  # coeffs = object containing coefficients and p-values of best model
  # title = title for graph (in quotes)
  # y_axis_label = label for y-axis (in quotes)
  ggplot(data, aes(x = date, y = {{ var }})) +
    geom_line(color = "burlywood") +
    geom_smooth(method = "lm", color = "black", size = 0.6) +
    labs(title = title,
         x = "Month",
         y = y_axis_label,
         caption = paste("slope = ", round(coeffs[2,1], 3), "\n(p-value = ", coeffs[2,2], ")"))
}


# functions to run nlme models ----------------------------------------------------------

# base model function
base_nlme_model <- function(data, var, month) {
  # data = dataframe
  # var = var to model
  # month = month of year
  
  data <- data %>% filter(month == month)
  
  model_specification = as.formula(paste0(var, " ~ time"))
  
  gls(model = model_specification,
      data = data,
      method = "ML",
      na.action = na.omit)
}



# AR1 model function
AR1_nlme_model <- function(data, var, month) {
  # data = dataframe
  # var = var to model
  # month = month of year
  
  data <- data %>% filter(month == month)
  
  model_specification = as.formula(paste0(var, " ~ time"))
  
  gls(model = model_specification,
      data = data,
      correlation = corARMA(form = ~ time, p = 1),
      method = "ML",
      na.action = na.omit)
}


# AR2 model function
AR2_nlme_model <- function(data, var, month) {
  # data = dataframe
  # var = var to model
  # month = month of year
  
  data <- data %>% filter(month == month)
  
  model_specification = as.formula(paste0(var, " ~ time"))
  
  gls(model = model_specification,
      data = data,
      correlation = corARMA(form = ~ time, p = 2),
      method = "ML",
      na.action = na.omit)
}

Mean Monthly Air Temperature -

plot_monthly_prelim(m, airt, "Monthly Mean Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

Met 40 - June:

by_month_lm(m, airt, "40", "Met 40 - Deep Well \nMean Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m40_6 <- prepare_data_for_nlme(m, "40", "6")

m40_6_mbase <- base_nlme_model(m40_6, "airt", 6)
m40_6_mAR1 <- AR1_nlme_model(m40_6, "airt", 6)
m40_5_mAR2 <- AR2_nlme_model(m40_6, "airt", 6)

AICc scores:

AICc(m40_6_mbase, m40_6_mAR1, m40_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m40_6_mbase 3 110.5176
m40_6_mAR1 4 112.6553
m40_5_mAR2 5 114.8095 
summary(m40_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   109.7176 114.2966 -51.85878
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 22.828769 0.4020514 56.78072  0.0000
## time         0.059215 0.0200400  2.95484  0.0058
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.8703743 -0.7313602  0.2242368  0.6785338  1.9307997 
## 
## Residual standard error: 1.112177 
## Degrees of freedom: 34 total; 32 residual
m40_6_coeffs <- summary(m40_6_mbase)$tTable[,c(1, 4)]
(m40_6_plot <- plot_monthly_results(m40_6, airt, m40_6_coeffs, "Met 40 - Deep Well \nJune Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 42 - June:

by_month_lm(m, airt, "42", "Met 42 - Deep Well \nMean Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m42_6 <- prepare_data_for_nlme(m, "42", "6")

m42_6_mbase <- base_nlme_model(m42_6, "airt", 6)
m42_6_mAR1 <- AR1_nlme_model(m42_6, "airt", 6)
m42_5_mAR2 <- AR2_nlme_model(m42_6, "airt", 6)

AICc scores:

AICc(m42_6_mbase, m42_6_mAR1, m42_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m42_6_mbase 3 105.9467
m42_6_mAR1 4 108.3905
m42_5_mAR2 5 111.1827 
summary(m42_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   105.1191 109.6086 -49.55955
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 21.896715 0.3992867 54.83958  0.0000
## time         0.046355 0.0204920  2.26211  0.0308
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.10520508 -0.68571065 -0.05950743  0.72953403  2.59403949 
## 
## Residual standard error: 1.086396 
## Degrees of freedom: 33 total; 31 residual
m42_6_coeffs <- summary(m42_6_mbase)$tTable[,c(1, 4)]
(m42_6_plot <- plot_monthly_results(m42_6, airt, m42_6_coeffs, "Met 42 - Cerro Montoso \nJune Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 49 - June:

by_month_lm(m49_6, airt, "49", "Met 49 - Deep Well \nMean Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m49_6 <- prepare_data_for_nlme(m, "49", "6")

m49_6_mbase <- base_nlme_model(m49_6, "airt", 6)
m49_6_mAR1 <- AR1_nlme_model(m49_6, "airt", 6)
m49_5_mAR2 <- AR2_nlme_model(m49_6, "airt", 6)

AICc scores:

AICc(m49_6_mbase, m49_6_mAR1, m49_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m49_6_mbase 3 71.53132
m49_6_mAR1 4 74.26520
m49_5_mAR2 5 77.55554 
summary(m49_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   70.26816 73.67465 -32.13408
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 23.986467 0.4413401 54.34917  0.0000
## time         0.066771 0.0321880  2.07440  0.0505
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2133523 -0.7132157  0.1294271  0.5771409  1.7943170 
## 
## Residual standard error: 0.9784315 
## Degrees of freedom: 23 total; 21 residual
m49_6_coeffs <- summary(m49_6_mbase)$tTable[,c(1, 4)]
(m49_6_plot <- plot_monthly_results(m49_6, airt, m49_6_coeffs, "Met 49 - Five Points \nJune Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 50 - June:

by_month_lm(m, airt, "50", "Met 50 - Deep Well \nMean Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m50_6 <- prepare_data_for_nlme(m, "50", "6")

m50_6_mbase <- base_nlme_model(m50_6, "airt", 6)
m50_6_mAR1 <- AR1_nlme_model(m50_6, "airt", 6)
m50_5_mAR2 <- AR2_nlme_model(m50_6, "airt", 6)

AICc scores:

AICc(m50_6_mbase, m50_6_mAR1, m50_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m50_6_mbase 3 68.93686
m50_6_mAR1 4 71.45522
m50_5_mAR2 5 75.03545 
summary(m50_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   67.43686 70.42406 -30.71843
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 25.066671 0.5504260 45.54049  0.0000
## time         0.072846 0.0459487  1.58538  0.1303
## 
##  Correlation: 
##      (Intr)
## time -0.877
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.30015742 -0.66646539  0.08760363  0.68926206  1.49495780 
## 
## Residual standard error: 1.1241 
## Degrees of freedom: 20 total; 18 residual
m50_6_coeffs <- summary(m50_6_mbase)$tTable[,c(1, 4)]
(m50_6_plot <- plot_monthly_results(m50_6, airt, m50_6_coeffs, "Met 50 - Blue Grama \nJune Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Mean Monthly Air Temperature for June

# airt results 
(airt_6_graphs <- grid.arrange(m40_6_plot, m42_6_plot, m49_6_plot, m50_6_plot, ncol=2))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]

Maximum Monthly Air Temperature -

plot_monthly_prelim(m, maxair, "Monthly Maximum Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

Met 40 - June:

by_month_lm(m, maxair, "40", "Met 40 - Deep Well \nMaximum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m40_6 <- prepare_data_for_nlme(m, "40", "6")

m40_6_mbase <- base_nlme_model(m40_6, "maxair", 6)
m40_6_mAR1 <- AR1_nlme_model(m40_6, "maxair", 6)
m40_5_mAR2 <- AR2_nlme_model(m40_6, "maxair", 6)

AICc scores:

AICc(m40_6_mbase, m40_6_mAR1, m40_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m40_6_mbase 3 131.0670
m40_6_mAR1 4 132.3726
m40_5_mAR2 5 134.7689 
summary(m40_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##       AIC      BIC    logLik
##   130.267 134.8461 -62.13351
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 36.74887 0.5439069 67.56463  0.0000
## time         0.07725 0.0271108  2.84947  0.0076
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.35174301 -0.43438926  0.07110273  0.48204954  2.65213862 
## 
## Residual standard error: 1.504586 
## Degrees of freedom: 34 total; 32 residual
m40_6_coeffs <- summary(m40_6_mbase)$tTable[,c(1, 4)]
(m40_6_plot <- plot_monthly_results(m40_6, maxair, m40_6_coeffs, "Met 40 - Deep Well \nJune Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 42 - June:

by_month_lm(m, maxair, "42", "Met 42 - Deep Well \nMaximum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m42_6 <- prepare_data_for_nlme(m, "42", "6")

m42_6_mbase <- base_nlme_model(m42_6, "maxair", 6)
m42_6_mAR1 <- AR1_nlme_model(m42_6, "maxair", 6)
m42_5_mAR2 <- AR2_nlme_model(m42_6, "maxair", 6)

AICc scores:

AICc(m42_6_mbase, m42_6_mAR1, m42_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m42_6_mbase 3 130.7457
m42_6_mAR1 4 133.1964
m42_5_mAR2 5 135.4923 
summary(m42_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   129.9181 134.4076 -61.95906
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 34.12709 0.5813903 58.69911  0.0000
## time         0.08010 0.0298378  2.68441  0.0116
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.76318154 -0.69378514 -0.02516334  0.53254899  3.12435430 
## 
## Residual standard error: 1.581872 
## Degrees of freedom: 33 total; 31 residual
m42_6_coeffs <- summary(m42_6_mbase)$tTable[,c(1, 4)]
(m42_6_plot <- plot_monthly_results(m42_6, maxair, m42_6_coeffs, "Met 42 - Cerro Montoso \nJune Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 49 - June:

by_month_lm(m49_6, maxair, "49", "Met 49 - Deep Well \nMaximum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m49_6 <- prepare_data_for_nlme(m, "49", "6")

m49_6_mbase <- base_nlme_model(m49_6, "maxair", 6)
m49_6_mAR1 <- AR1_nlme_model(m49_6, "maxair", 6)
m49_5_mAR2 <- AR2_nlme_model(m49_6, "maxair", 6)

AICc scores:

AICc(m49_6_mbase, m49_6_mAR1, m49_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m49_6_mbase 3 82.12418
m49_6_mAR1 4 84.65632
m49_5_mAR2 5 87.96102 
summary(m49_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   80.86102 84.26751 -37.43051
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 37.52901 0.5556259 67.54367  0.0000
## time         0.05914 0.0405232  1.45942  0.1592
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.5406398 -0.3613768  0.2014793  0.6971315  1.8703417 
## 
## Residual standard error: 1.231798 
## Degrees of freedom: 23 total; 21 residual
m49_6_coeffs <- summary(m49_6_mbase)$tTable[,c(1, 4)]
(m49_6_plot <- plot_monthly_results(m49_6, maxair, m49_6_coeffs, "Met 49 - Five Points \nJune Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 50 - June:

by_month_lm(m, maxair, "50", "Met 50 - Deep Well \nMaximum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m50_6 <- prepare_data_for_nlme(m, "50", "6")

m50_6_mbase <- base_nlme_model(m50_6, "maxair", 6)
m50_6_mAR1 <- AR1_nlme_model(m50_6, "maxair", 6)
m50_5_mAR2 <- AR2_nlme_model(m50_6, "maxair", 6)

AICc scores:

AICc(m50_6_mbase, m50_6_mAR1, m50_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m50_6_mbase 3 78.64979
m50_6_mAR1 4 81.79530
m50_5_mAR2 5 85.41433 
summary(m50_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   77.14979 80.13699 -35.57489
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 37.54153 0.7017069 53.50030  0.0000
## time         0.12119 0.0585774  2.06885  0.0532
## 
##  Correlation: 
##      (Intr)
## time -0.877
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.3244313 -0.5513157  0.1981575  0.5670821  1.9498376 
## 
## Residual standard error: 1.433052 
## Degrees of freedom: 20 total; 18 residual
m50_6_coeffs <- summary(m50_6_mbase)$tTable[,c(1, 4)]
(m50_6_plot <- plot_monthly_results(m50_6, maxair, m50_6_coeffs, "Met 50 - Blue Grama \nJune Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Maximum Monthly Air Temperature for June

# maxair results 
(maxair_6_graphs <- grid.arrange(m40_6_plot, m42_6_plot, m49_6_plot, m50_6_plot, ncol=2))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]

Minimum Monthly Air Temperature -

plot_monthly_prelim(m, minair, "Monthly Minimum Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

Met 40 - June:

by_month_lm(m, minair, "40", "Met 40 - Deep Well \nMinimum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m40_6 <- prepare_data_for_nlme(m, "40", "6")

m40_6_mbase <- base_nlme_model(m40_6, "minair", 6)
m40_6_mAR1 <- AR1_nlme_model(m40_6, "minair", 6)
m40_5_mAR2 <- AR2_nlme_model(m40_6, "minair", 6)

AICc scores:

AICc(m40_6_mbase, m40_6_mAR1, m40_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m40_6_mbase 3 168.6741
m40_6_mAR1 4 171.2504
m40_5_mAR2 5 173.8996 
summary(m40_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   167.8741 172.4532 -80.93704
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept)  7.833538 0.9456032  8.284170  0.0000
## time        -0.029816 0.0471331 -0.632596  0.5315
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.25544260 -0.73370381 -0.09415726  0.75553571  1.62187112 
## 
## Residual standard error: 2.615781 
## Degrees of freedom: 34 total; 32 residual
m40_6_coeffs <- summary(m40_6_mbase)$tTable[,c(1, 4)]
(m40_6_plot <- plot_monthly_results(m40_6, minair, m40_6_coeffs, "Met 40 - Deep Well \nJune Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 42 - June:

by_month_lm(m, minair, "42", "Met 42 - Deep Well \nMinimum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m42_6 <- prepare_data_for_nlme(m, "42", "6")

m42_6_mbase <- base_nlme_model(m42_6, "minair", 6)
m42_6_mAR1 <- AR1_nlme_model(m42_6, "minair", 6)
m42_5_mAR2 <- AR2_nlme_model(m42_6, "minair", 6)

AICc scores:

AICc(m42_6_mbase, m42_6_mAR1, m42_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m42_6_mbase 3 153.6229
m42_6_mAR1 4 155.3432
m42_5_mAR2 5 155.8929 
summary(m42_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   152.7953 157.2848 -73.39764
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept)  8.669623 0.8222512 10.54376  0.0000
## time        -0.039671 0.0421991 -0.94009  0.3544
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2567953 -0.6711610 -0.1724796  0.7295458  2.1228551 
## 
## Residual standard error: 2.237216 
## Degrees of freedom: 33 total; 31 residual
m42_6_coeffs <- summary(m42_6_mbase)$tTable[,c(1, 4)]
(m42_6_plot <- plot_monthly_results(m42_6, minair, m42_6_coeffs, "Met 42 - Cerro Montoso \nJune Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 49 - June:

by_month_lm(m49_6, minair, "49", "Met 49 - Deep Well \nMinimum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m49_6 <- prepare_data_for_nlme(m, "49", "6")

m49_6_mbase <- base_nlme_model(m49_6, "minair", 6)
m49_6_mAR1 <- AR1_nlme_model(m49_6, "minair", 6)
m49_5_mAR2 <- AR2_nlme_model(m49_6, "minair", 6)

AICc scores:

AICc(m49_6_mbase, m49_6_mAR1, m49_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m49_6_mbase 3 103.5652
m49_6_mAR1 4 101.6278
m49_5_mAR2 5 104.2407 
summary(m49_6_mAR1)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   99.40558 103.9476 -45.70279
## 
## Correlation Structure: AR(1)
##  Formula: ~time 
##  Parameter estimate(s):
##        Phi 
## -0.4759774 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 8.610738 0.5561496 15.482773  0.0000
## time        0.058075 0.0408268  1.422464  0.1696
## 
##  Correlation: 
##      (Intr)
## time -0.881
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.97613576 -0.57321163 -0.02887857  0.77814343  1.67707143 
## 
## Residual standard error: 1.995719 
## Degrees of freedom: 23 total; 21 residual
m49_6_coeffs <- summary(m49_6_mAR1)$tTable[,c(1, 4)]
(m49_6_plot <- plot_monthly_results(m49_6, minair, m49_6_coeffs, "Met 49 - Five Points \nJune Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 50 - June:

by_month_lm(m, minair, "50", "Met 50 - Deep Well \nMinimum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

m50_6 <- prepare_data_for_nlme(m, "50", "6")

m50_6_mbase <- base_nlme_model(m50_6, "minair", 6)
m50_6_mAR1 <- AR1_nlme_model(m50_6, "minair", 6)
m50_5_mAR2 <- AR2_nlme_model(m50_6, "minair", 6)

AICc scores:

AICc(m50_6_mbase, m50_6_mAR1, m50_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m50_6_mbase 3 90.80023
m50_6_mAR1 4 90.42714
m50_5_mAR2 5 93.73701 
summary(m50_6_mAR1)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   87.76048 91.74341 -39.88024
## 
## Correlation Structure: AR(1)
##  Formula: ~time 
##  Parameter estimate(s):
##       Phi 
## 0.4078502 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 12.354002 1.3459652  9.178545   0.000
## time        -0.104734 0.1105605 -0.947298   0.356
## 
##  Correlation: 
##      (Intr)
## time -0.862
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.9317235 -0.5736404  0.1539656  0.5339572  2.3574953 
## 
## Residual standard error: 1.937679 
## Degrees of freedom: 20 total; 18 residual
m50_6_coeffs <- summary(m50_6_mAR1)$tTable[,c(1, 4)]
(m50_6_plot <- plot_monthly_results(m50_6, minair, m50_6_coeffs, "Met 50 - Blue Grama \nJune Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Minimum Monthly Air Temperature for June

# minair results 
(minair_6_graphs <- grid.arrange(m40_6_plot, m42_6_plot, m49_6_plot, m50_6_plot, ncol=2))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]

June Minimum, Mean, and Maximum Air Temperature Results -

minair_6_graphs
## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
airt_6_graphs
## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
maxair_6_graphs
## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]

Monthly Precipitation -

plot_monthly_prelim(m, ppt, "Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'

Met 40 - June:

by_month_lm(m, ppt, "40", "Met 40 - Deep Well \n Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'

m40_6 <- prepare_data_for_nlme(m, "40", "6")

m40_6_mbase <- base_nlme_model(m40_6, "ppt", 6)
m40_6_mAR1 <- AR1_nlme_model(m40_6, "ppt", 6)
m40_5_mAR2 <- AR2_nlme_model(m40_6, "ppt", 6)

AICc scores:

AICc(m40_6_mbase, m40_6_mAR1, m40_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m40_6_mbase 3 303.5660
m40_6_mAR1 4 306.0801
m40_5_mAR2 5 306.6254 
summary(m40_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##       AIC      BIC   logLik
##   302.766 307.3451 -148.383
## 
## Coefficients:
##                 Value Std.Error    t-value p-value
## (Intercept) 21.708889  6.874179  3.1580337  0.0035
## time        -0.282118  0.342640 -0.8233652  0.4164
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.1119553 -0.7733032 -0.3014592  0.3845706  2.9864295 
## 
## Residual standard error: 19.01574 
## Degrees of freedom: 34 total; 32 residual
m40_6_coeffs <- summary(m40_6_mbase)$tTable[,c(1, 4)]
(m40_6_plot <- plot_monthly_results(m40_6, ppt, m40_6_coeffs, "Met 40 - Deep Well \nJune Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'

Met 42 - June:

by_month_lm(m, ppt, "42", "Met 42 - Deep Well \n Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'

m42_6 <- prepare_data_for_nlme(m, "42", "6")

m42_6_mbase <- base_nlme_model(m42_6, "ppt", 6)
m42_6_mAR1 <- AR1_nlme_model(m42_6, "ppt", 6)
m42_5_mAR2 <- AR2_nlme_model(m42_6, "ppt", 6)

AICc scores:

AICc(m42_6_mbase, m42_6_mAR1, m42_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m42_6_mbase 3 303.7262
m42_6_mAR1 4 306.0666
m42_5_mAR2 5 305.4861 
summary(m42_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   302.8986 307.3881 -148.4493
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 26.517472  7.993137  3.317530  0.0023
## time        -0.355085  0.410219 -0.865599  0.3934
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.2029751 -0.6962164 -0.3386889  0.3329534  3.9586599 
## 
## Residual standard error: 21.74807 
## Degrees of freedom: 33 total; 31 residual
m42_6_coeffs <- summary(m42_6_mbase)$tTable[,c(1, 4)]
(m42_6_plot <- plot_monthly_results(m42_6, ppt, m42_6_coeffs, "Met 42 - Cerro Montoso \nJune Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'

Met 49 - June:

by_month_lm(m49_6, ppt, "49", "Met 49 - Deep Well \n Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'

m49_6 <- prepare_data_for_nlme(m, "49", "6")

m49_6_mbase <- base_nlme_model(m49_6, "ppt", 6)
m49_6_mAR1 <- AR1_nlme_model(m49_6, "ppt", 6)
m49_5_mAR2 <- AR2_nlme_model(m49_6, "ppt", 6)

AICc scores:

AICc(m49_6_mbase, m49_6_mAR1, m49_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m49_6_mbase 3 177.3561
m49_6_mAR1 4 178.3878
m49_5_mAR2 5 175.5247 
summary(m49_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   176.0929 179.4994 -85.04645
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 11.236151  4.404374 2.5511347  0.0186
## time         0.111922  0.321222 0.3484245  0.7310
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.2997479 -0.8456391 -0.3466543  0.7219178  2.1507099 
## 
## Residual standard error: 9.764303 
## Degrees of freedom: 23 total; 21 residual
m49_6_coeffs <- summary(m49_6_mbase)$tTable[,c(1, 4)]
(m49_6_plot <- plot_monthly_results(m49_6, ppt, m49_6_coeffs, "Met 49 - Five Points \nJune Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'

Met 50 - June:

by_month_lm(m, ppt, "50", "Met 50 - Deep Well \n Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'

m50_6 <- prepare_data_for_nlme(m, "50", "6")

m50_6_mbase <- base_nlme_model(m50_6, "ppt", 6)
m50_6_mAR1 <- AR1_nlme_model(m50_6, "ppt", 6)
m50_5_mAR2 <- AR2_nlme_model(m50_6, "ppt", 6)

AICc scores:

AICc(m50_6_mbase, m50_6_mAR1, m50_5_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m50_6_mbase 3 176.3212
m50_6_mAR1 4 177.2593
m50_5_mAR2 5 178.6292 
summary(m50_6_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   174.8212 177.8084 -84.41058
## 
## Coefficients:
##                 Value Std.Error    t-value p-value
## (Intercept) 21.677368  8.065088  2.6878030  0.0150
## time        -0.611654  0.673261 -0.9084954  0.3756
## 
##  Correlation: 
##      (Intr)
## time -0.877
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.0704234 -0.5841853 -0.2183698  0.4232182  3.4459480 
## 
## Residual standard error: 16.47082 
## Degrees of freedom: 20 total; 18 residual
m50_6_coeffs <- summary(m50_6_mbase)$tTable[,c(1, 4)]
(m50_6_plot <- plot_monthly_results(m50_6, ppt, m50_6_coeffs, "Met 50 - Blue Grama \nJune Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Monthly Precipitation for June

# ppt results 
(ppt_6_graphs <- grid.arrange(m40_6_plot, m42_6_plot, m49_6_plot, m50_6_plot, ncol=2))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]